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Abstract 
CN . 

^ 'This work addresses computing techniques for dose calculations in treatment planning with proton and ion beams, based 
^ 'on an efficient kernel-convolution method referred to as grid-dose spreading (GDS) and accurate heterogeneity-correction 

method referred to as Gaussian beam splitting. The original GDS algorithm suffered from distortion of dose distribution 
\0 !for beams tilted with respect to the dose-grid axes. Use of intermediate grids normal to the beam field has solved 

'the beam-tilting distortion. Interplay of arrangement between beams and grids was found as another intrinsic source 

'of artifact. Inclusion of rectangular- kernel convolution in beam transport, to share the beam contribution among the 
f~| 'nearest grids in a regulatory manner, has solved the interplay problem. This algorithmic framework was applied to a 
On'tilted proton pencil beam and a broad carbon-ion beam. In these cases, while the elementary pencil beams individually 

split into several tens, the calculation time increased only by several times with the GDS algorithm. The GDS and bcam- 
^ splitting methods will complementarily enable accurate and efficient dose calculations for radiotherapy with protons and 
^ 'ions. 

ty^ •Keywords: proton, ion beam, pencil-beam algorithm, treatment planning, inhomogeneity correction 
_0 ]PACS: 87.55.D-, 87.55.kd 
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.1. Introduction 

Dose distributions of radiotherapy are represented by 
point doses at orthogonally arranged grids. In treatment- 
planning practice, the grid intervals are defined from a 
physical, clinical, and practical points of view, often re- 
sulting in cubic dimensions of a few millimeters. Accuracy, 
efficiency and their balance are essential in practice, for 
which the pencil-beam algorithm is commonly used. That 
is mathematically a convolution integral of total energy 
released per mass (terma) with elementary beam-spread 
ikernel, which may be computationally demanding. 

The grid-dose-spreading (GDS) algorithm was devel- 
pped for fast dose calculation of heavy-charged-particle 
■beams in patient body [1]. The GDS algorithm employs 
approximation to extract beam-interaction part from the 
■integral at the expense of distortion of dose distribution 
for a beam tilted with respect to the grid axes, as origi- 
nally recognized in Ref. [1]. The beam-tilting distortion 
may be generally insignificant when beam blurring is as 
small as the required spatial resolution, for example, for 
a carbon-ion beam. In fact, the GDS method was suc- 
cessfully incorporated into a clinical treatment-planning 
system for carbon-ion radiotherapy with vertical and hori- 
zontal fixed beams [2, 3], for which tilting was intrinsically 
absent. 
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In that particular implementation, a simplistic post 
process was added to the original broad-beam algorithm 
so as to spread an intermediate terma distribution uni- 
formly [1]. In general, the spreading kernel could be spa- 
tially modulated using the pencil-beam model for more 
accurate heterogeneity correction [4]. There are two re- 
ciprocal approaches for convolution, i.e. to collect doses 
transferred from nearby interactions to a grid or the dose- 
deposition point of view and to spread a terma from an in- 
teraction to nearby grids or the interaction point of view. 
The latter is usually more efficient than the former for 
three-dimensional dose calculation [5]. 

The pencil-beam model implicitly assumes homogene- 
ity of the medium within the elementary beam spread. 
Beams that have grown excessively thick in heterogeneous 
transport are thus incompatible. As a general and rigor- 
ous solution, Gaussian-beam splitting was proposed, with 
which overgrown beams are subdivided into smaller ones 
at locations of large lateral heterogeneity [6]. Figure 1 
demonstrates its effectiveness for a simple density bound- 
ary, where the non-splitting beam happened to traverse 
an edge of a bone-equivalent material while about a half of 
the split beams traverse the bone-equivalent material. The 
splitting causes explosive beam multiplication in a shower- 
like process. In this particular case for example, the orig- 
inal beam recursively split into 28 final beams. Slowing 
down of dose calculation due to beam multiplication will 
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Figure 1: (a) Non-splitting and (b) splitting dose calculations with 
isodose lines at every 10% levels of the maximum non-splitting dose 
in the y = cross section, where a proton pencil beam with E = 150 
MeV and cr = 3 mm is incident into water with a bone-equivalent 
material {p = 1.8 g/cm'^) inserted halfway (gray area). 

be a problem in practice. 

In Ref. [6], the beam-splitting method was stated as 
efficient due to certain "algorithmic techniques to be ex- 
plained elsewhere" , which in fact implied this work to con- 
struct a framework, where the GDS and beam-splitting 
methods work compatibly for accurate and efficient dose 
calculations. In addition, we will refine the GDS algorithm 
with a fix against the beam-tilting distortion and with the 
pencil-beam model in the interaction point of view for bet- 
ter heterogeneity correction. 

Although the Gaussian-beam approximation may be 
reasonable for the multiple-scattering effect, two or more 
Gaussian components would improve the accuracy of lat- 
eral dose distribution of proton and ion pencil beams [7, 8] . 
However, such large-sized components are intrinsically in- 
compatible with fine heterogeneity. In addition, it is incon- 
ceivable to apply the beam-splitting method for large-sized 
components to secure practical efficiency. 

This framework will be applicable not only to broad- 
beam delivery but also to pencil-beam scanning, where a 
physical scanned beam may have to be decomposed into 
virtual elementary beams to address heterogeneity [9] . As 
this work aims to improve computing methods, we focus 
on evaluation of efficiency and settlement of the intrinsic 
artifacts with respect to the ideal beam models that are 
mathematically given, without repeating experimental as- 
sessments of accuracy [6] . 

2. Materials and methods 

2.1. Algorithmic techniques 
2.1.1. Grid normalization 

We will solve the beam-tilting distortion of the GDS 
algorithm by defining intermediate grids for dose calcula- 
tion, which are arranged to be normal to the beam-field 
axes. As shown in Figure 2, the original dose grids along 



Figure 2: Schematic of grid normalization in cross-section view per- 
pendicular to 63 = , where shown are a field (gray area) with axes 
Sx and Cz, margins (light gray areas), original grids with axes ei and 
62, and normal grids of width W and height H. 

numbered axes 1,2, and 3 are defined with basis vectors 
ei, 6*2, and f?3 and intervals 82, and (S3. For a given 
radiation field, the field coordinates x, y, and z with basis 
vectors e^, Cy, and Cz are associated, where the origin is at 
the isocenter and Cz is in the source direction. With lateral 
margins for penumbra, the normal-grid volume is defined 
as the suprcmum of normal rectangular-parallelepiped vol- 
ume oi W X L X H containing the original grids in the 
margined field. Quadratic projection of the original-grid 
voxel gives the normal-grid intervals Sx, 6y, and Sz as 
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^h,y,^} = ^^a (ea ■ e{x,y,z}) , (1) 

a=l 

to approximately conserve the equivalent resolution. Nor- 
mal grids gijk are defined at equally spaced positions fg. 
for indices i € [0, [W'^rr]], j G [0, [L/Sy]] and k G [0, \hJ5z]], 
where [ ] is the ceiling function. 

Because the elementary pencil beams are almost par- 
allel to vector ez , normal-grid doses Dg can be accurately 
and efficiently calculated with 

rys+Sy/2 

^9 = Y.^h G,-,,,- (a;)da: / Gj,- (y) dy, 

^ JXg-5^/2 Jyg-Sy/2 

(2) 

where T^, cr^, and {xf^,yf^) are the terma, the spread, and 
the position of normal grid h, and Grn,a{x) is the Gaussian 
distribution of mean m and standard deviation a whose 
integral is readily given by the standard error function. 
Dose Dg at original grid g is then given by trilinear inter- 
polation, which is repeated for all the original grids in the 
margined field. In this manner, we only deal with normal 
grids in the GDS algorithm hereafter. 

2.1.2. Inter- grid beam sharing 

Beam-6 spread ab,g at grid g is given by beam-transport 
calculation [10, 11] and terma contribution ATi,,g is defined 
[6] as 



where Ub is the number of particles that is modeled as a 
constant, I)$o is the dose per in-air fluence that is mea- 
sured as a depth-dose curve, and As is the step length in 
voxel g. In the original GDS algorithm, grid terma Tg and 
spread ag are intermediately defined with formulas 
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(4) 



which arc straightforward extensions of Eqs. (13) and (14) 
in Ref. [1]. That would be, however, problematic when the 
beams are arranged independently to the grids. As shown 
in Figure 3(a), only several beams may traverse each voxel. 
Due to interplay between the beam and grid structures, 
the total beam-path length per voxel and consequently the 
grid terma Tg largely fluctuates. The fluctuation would be 
smeared out only with large spreading with ag ^ ^{x,y}- 

To fully resolve the fluctuation, we will distribute step 
terma ATf,^g to the nearest four grids in a regulatory man- 
ner. As shown in Figures 3(b) and 3(c), the x and y axes 
comprise the lateral plane, where any beam b is now mod- 
eled to have cross section Ah of size 6x x Sy centered at 
mid-step point fi, = {xt, yt, z&) in current voxel V/j at grid 
h. The four voxels that intersect Af, share terma contribu- 
tion ATb^h by areal fractions. Equation (4) is then modi- 
fied in such a way that the terma and spread distributions 
will be formed when all the beams have been processed as 
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where the areal fraction is given by 
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(5) 
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(7) 



for nearby grid g with |x6 — a;g| < 5x and \yb — yg\ < or 
otherwise 0. 

Essentially, this beam-sharing operation is one form 
of convolution with a rectangular kernel. Wc thus modify 
axial dose spreading a^, at grid h in (2) to correct the extra 
rectangular spreading as 
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(8) 



for which, sufficiently fine grid intervals must be 

given. 

2.2. Evaluation 

Tilted pencil beam. We examine the effect of grid normal- 
ization using an analytic 150 McV proton pencil beam 
model with formulated dose per in- air fluence [12] and 
lateral spread [11]. The beam with rms projected angle 
00 = 10 mrad and size (Tq = 1-2 mm was placed at the 
isocenter to incident into water in the direction of patient- 
support angle 6s = 0° and gantry angle (/)g = 30° in the 




Figure 3: Schematic of inter-grid beam sharing in the cross-section 
views with grid points (+), voxel boundaries (dashed hnes), beam 
paths (arrows), mid-step points (®), and rectangular A5 (gray ar- 
eas). One of many beams in (a) is focused in (b), where the dotted 
line indicates the x-y cross section in (c). 



standard coordinate system [13]. The original grids were 
defined cubically as (5i ~ (52 = (^3 = 1 mm to form a 
volume of 10 cm x 15 cm x 10 cm. Around the beam 
axis, 1-cm margins were added to define the normal-grid 
volume. The dose distributions were calculated using the 
GDS algorithm with and without grid normalization for 
comparison. 

Broad beam with heterogeneity. We examine a broad beam 
in a heterogeneous medium, which approximates clinical 
situations. The dose per in-air fluence was modeled based 
on an experiment [6], where a carbon-ion beam with inci- 
dent nucleon kinetic energy E/A = 290 MeV was broad- 
ened by spiral wobbling [14] and was range-modulated by 
a semi-Gaussian filter of an = 1.8 mm. In the calculation, 
a radiation from a source at height 500 cm was limited by 
a collimator at height 50 cm to form a 4 cm x 4 cm field on 
the isocenter plane, which yielded 81 x 81 original pencil 
beams. In a 20 cm x 20 cm x 10 cm volume griddcd at in- 
tervals of 1 mm along axes 1,2, and 3, an 19 cm diameter 
cylindrical water phantom was defined at the center. The 
phantom included a 2 cm diameter dense {p = 2 g/cm^) 
rod at radius 6 cm. The dosimetric effects of grid normal- 
ization and inter-grid sharing were investigated for two 
gantry angles 4>g = 0° and 45° with the concurrently ro- 
tated phantom so that the rod was always in the middle of 
the field to see any angular artifact. The computing times 
were measured using 2.4 GHz Intel Core2Duo processor on 
Apple MacBook computer. 

3. Results 

Tilted pencil beam. Figure 4 shows the projected dose dis- 
tributions for the tilted proton pencil beam. The original 
GDS algorithm severely distorted the analytic model that 
was exactly the dose kernel of the algorithm. The grid 
normalization greatly reduced the distortion. 

Heterogeneous broad beam. Table 1 and Figures 5 and 6 
show the configurations for broad carbon-ion beam calcu- 
lations, the computing times, and the dose distribution in 
the y = cross-section and along the x and z axes of the 
beam field. Configurations A and D resulted in severe dose 
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Figure 4: Projected dose at every 10% levels for (a) the analytic 
proton-beam model, (b) the original GDS calculation, and (c) the 
normalized GDS calculation. 



Table 1: Definitions of calculation configurations A-F in combina- 
tions of with (-f) or without {— ) 45° tilt T, grid normalization A'', 
intcr-grid sharing /, and beam splitting S and total computing times 
for a broad beam in a heterogeneous phantom. 



Config. T N I S Time (s) 



A 
B 
C 
D 
E 
F 



1.9 
2.0 
15.5 
3.0 
2.5 
16.4 
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Figure 5: Dose distributions in the y = cross section for a carbon- 
ion beam calculated in configurations A-F. The solid and dotted 
lines indicate every 10% dose levels and the field x and z axes. The 
light-gray and medium-gray areas indicate the water phantom and 
the high-density rod. 



artifacts for the beam-grid interplay and for the tilted in- 
cidence while the others are hardly distinguishable. From 
the A-B comparison, we find that the inter-grid sharing 
resolved the interplay artifact and added marginal com- 
puting load. From the D-E comparison, we find that the 
grid normalization resolved the angle issue and unexpect- 
edly reduced the computing time. From the B-C and E- 
F comparisons, we find that the beam splitting increased 
the computing time by factors 7.8 and 6.6 although the 
dosimetric effect was only marginal dose fluctuation from 
moderate heterogeneity of the round structures. 

4. Discussion 

In this work, wc introduced grid normalization to suc- 
cessfully resolve the problems with tilted incidence. The 
interpolation errors are limited by the grid intervals, which 
may be generally tolerable and controllable because the 
grid intervals are normally specified by a treatment plan- 
ner. We found no apparent loss of speed for interpolation. 
In fact, the areal convolution with normal grids was faster 
than the volumetric convolution with tilted grids. 

The beam-grid interplay artifact was overlooked in the 
original formulation of the GDS algorithm [1] that was 
implemented as an extension to the broad-beam algorithm 
so that the termas and the spreads were calculated exactly 
at the grids and were free from the interplay. The inter- 
grid beam sharing, which is essentially rectangular-kernel 
convolution, has fully resolved the problem. 
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Figure 6: Dose profiles along the (a) x and (b) z axes for a carbon- 
ion beam in configurations A (solid), B (dashed), C (dotted), D (-I-), 
E (o), and F (x) in units of arbitrary reference dose U. 
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As shown in Ref. [6], the beam sphtting resolved the 
intrinsic problems with the pencil-beam model for fine het- 
erogeneity by dynamic subdivision. This inevitably re- 
quires the interaction-point-of-view approach because oth- 
erwise it would be difficult to trace histories of the split 
beams backwards from each dose-deposition grid [5] . Since 
the GDS convolution is directly coupled not to individual 
beams but to resultant tcrma and spread distributions, 
beam multiplication due to splitting will not increase the 
computing time for dose convolution. In other words, the 
combination of the GDS and beam-splitting methods is a 
rational consequence. 

The observed slowing by a factor of several times was 
mainly attributed to increased transport calculation of the 
split beams. These slowing factors were comparable to 
that originally reported [6], not surprisingly because they 
in fact used a similar framework of the beam-splitting GDS 
algorithm though without grid normalization nor inter- 
grid sharing. Those examples happened to be of normal 
incidence and of small interplay and were only intended 
for a proof of principle of beam splitting. 

The tcrma-wcightcd mean for the gridded beam spread 
(jg in Eq. (6) can not be generally valid because the mean 
spread only approximately represents the contributing Gaus- 
sian components. It was originally assumed that its vari- 
ation would be small enough to be handled as locally uni- 
form. However, beam splitting changes the spread abruptly. 
On the other hand, the beam splitting converts the ma- 
jor part of beam spreads directly into a terma distribution. 
As a result, grid-dose spreading handles only the residuals, 
for which the terma-weighted mean may suffice. 

5. Conclusions 

A known problem of tilted incidence with the origi- 
nal GDS algorithm was naturally resolved by the grid- 
normalization method without serious loss of accuracy or 
efficiency. Another problem of beam-grid interplay arti- 
fact was revealed and was resolved by the inter-grid beam- 
sharing method. 

The beam-splitting method for fine-heterogeneity cor- 
rection will inevitably multiply beams to transport and 
thus will slow down dose calculation. With the GDS al- 
gorithm, the dose convolution is made only once after all 
the beams have been transported, which minimizes the 
impact of the beam multiplication on computing time. In 
fact, for the beams individually split into several tens, the 
calculation time increased only by several times with the 
GDS. This algorithmic framework will thus enable fast and 
accurate treatment planning of heavy charged particle ra- 
diotherapy in the presence of density heterogeneity finer 
than the size of intrinsic beam blurring. 
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